home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / r_test.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  117 lines
  1. ;$Id: r_test.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       R_TEST
  8. ;
  9. ; PURPOSE:
  10. ;       This function tests the hypothesis that a binary sequence (a 
  11. ;       sequence of 1s and 0s) represents a "random sampling". This
  12. ;       test is based on the "theory of runs" and is often refered to
  13. ;       as the Runs Test for Randomness. The result is a two-element
  14. ;       vector containing the nearly-normal test statistic Z and the 
  15. ;       one-tailed probability of obtaining a value of Z or greater.
  16. ;
  17. ; CATEGORY:
  18. ;       Statistics.
  19. ;
  20. ; CALLING SEQUENCE:
  21. ;       Result = R_test(X)
  22. ;
  23. ; INPUTS:
  24. ;       X:    An n-element vector of type integer, float or double.
  25. ;             Elements not equal to 0 or 1 are removed and the length
  26. ;             of X is correspondingly reduced.
  27. ;
  28. ; KEYWORD PARAMETERS:
  29. ;       R:    Use this keyword to specify a named variable which returns
  30. ;             the number of runs (clusters of 0s and 1s) in X.       
  31. ;
  32. ;      N0:    Use this keyword to specify a named variable which returns
  33. ;             the number of 0s in X.
  34. ;
  35. ;      N1:    Use this keyword to specify a named variable which returns
  36. ;             the number of 1s in X.
  37. ;
  38. ; EXAMPLE:
  39. ;       Define a vector of 1s and 0s.
  40. ;         x = [0,1,1,0,1,0,0,0,1,0,0,1,1,0,1,0,1,0,0,1,0,1,1,0,1,0,0,1,0,1]
  41. ;
  42. ;       Test the hypothesis that x represents a random sampling against the 
  43. ;       hypothesis that it does not represent a random sampling at the 0.05 
  44. ;       significance level.
  45. ;         result = r_test(x, r = r, n0 = n0, n1 = n1)
  46. ;
  47. ;       The result should be the 2-element vector:
  48. ;         [2.26487, 0.0117604]
  49. ;       The keyword parameters should be returned as:
  50. ;         r = 22.0000, n0 = 16.0000, n1 = 14.0000 
  51. ;
  52. ;       The computed probability (0.0117604) is less than the 0.05
  53. ;       significance level and therefore we reject the hypothesis that x
  54. ;       represents a random sampling. The results show that there are too 
  55. ;       many runs, indicating a non-random cyclical pattern.
  56. ;
  57. ; PROCEDURE:
  58. ;       R_TEST computes the nonparametric Runs Test for Randomness. A 
  59. ;       "run" is a cluster of identical symbols within a sequence of two
  60. ;       distinct symbols. The binary sequence (x) defined in EXAMPLE has 
  61. ;       22 runs (or clusters). The first run contains one 0, the second 
  62. ;       run contains two 1s, the third run contains one 0, and so on.
  63. ;       In general, the randomness hypothesis will be rejected if there
  64. ;       are lengthy runs of 0s or 1s or if there are alternating patters
  65. ;       of many short runs.
  66. ;
  67. ; REFERENCE:
  68. ;       PROBABILITY and STATISTICS for ENGINEERS and SCIENTISTS (3rd edition)
  69. ;       Ronald E. Walpole & Raymond H. Myers
  70. ;       ISBN 0-02-424170-9
  71. ;
  72. ; MODIFICATION HISTORY:
  73. ;       Written by:  GGS, RSI, August 1994
  74. ;-
  75.  
  76. function r_test, x, r = r, n0 = n0, n1 = n1
  77.  
  78.   on_error, 2
  79.  
  80.   nx = n_elements(x)
  81.  
  82.   if nx le 10 then message, $
  83.     'Not defined for input vectors of 10 or fewer elements.'
  84.  
  85.   ;Remove any sequence elements that are not 0s or 1s.
  86.   data = where(x eq 0 or x eq 1, nb)
  87.   if nb ne 0 then x = x[data]
  88.  
  89.   sx = size(x)
  90.  
  91.   h0 = where(x eq 0, n0)
  92.   if n0 ne 0 then hi = where(h0+1 ne shift(h0, -1), nn0) $
  93.   else nn0 = 0L
  94.  
  95.   n1 = sx[1] - n0
  96.  
  97.   if x[sx[1]-1] ne x[0] then nn1 = nn0 $
  98.   else if x[0] eq 1 then nn1 = nn0 + 1 $
  99.   else nn1 = nn0 - 1
  100.  
  101.   if n0 eq 0 or n1 eq 0 then message, $
  102.     'x is a sequence of identical data.'
  103.  
  104.   r = nn0 + nn1
  105.   e = 2.0 * n0 * n1 / (n0 + n1) + 1.0
  106.   v = 2.0 * n0 * n1 * (2.0 * n1 * n0 - n0 - n1) / $
  107.                       ((n0 + n1 -1.0) * (n1 + n0)^2)
  108.   z = (r - e) / sqrt(v)
  109.  
  110.   prob = 1.0 - gauss_pdf(abs(z))
  111.  
  112.   return, [z, prob]
  113.  
  114. end
  115.  
  116.  
  117.